#ifndef __REFLOAT__
#define __REFLOAT__

#include <stdio.h>
#include <stdlib.h>
#include <mpfr.h>
#include <string>
#include <vector>
#include <iostream>
#include <stdexcept>

typedef unsigned long int ulong;

class refloat
{
  public:
    refloat(std::nullptr_t ptr);
    refloat(){}
    refloat(long int re, mpfr_prec_t prec=4, bool temporary=false);
    refloat(std::string re, mpfr_prec_t prec=4, int base=16, bool temporary=false);
    refloat(const refloat& z);
    ~refloat();
    void set_temporary(mpfr_prec_t prec);
    void set(const refloat& x);
    void set(long int re, mpfr_prec_t prec);
    void set(std::string re, mpfr_prec_t prec, int base);
    void set(long int re);
    void set(std::string re, int base);
    void set_prec(mpfr_prec_t prec);
    mpfr_prec_t get_prec() const;
    void neg();
    void add(const refloat& x);
    void add(long int i);
    void sub(const refloat& x);
    void sub(long int i);
    void mul(const refloat& x);
    void sqr();
    void div(const refloat& x);
    void sqrt();
    void mul_power_2(long int p);
    void reciprocal();
    void round();
    void abs();
    void print(std::string &re, int base=16, mpfr_prec_t prec=0, bool bold=true) const;
    bool isTemp() const;
    bool isSet() const;
    void prefetch() const;
    unsigned long int to_ulong();
    
    bool is_nan();
    
    refloat& operator =(const refloat& b);
    refloat operator +(const refloat& b) const;
    refloat operator -(const refloat& b) const;
    refloat operator *(const refloat& b) const;
    refloat operator /(const refloat& b) const;
    bool operator ==(const refloat& b) const;
    bool operator !=(const refloat& b) const;
    bool operator >=(const refloat& b) const;
    bool operator <=(const refloat& b) const;
    bool operator >(const refloat& b) const;
    bool operator <(const refloat& b) const;
    refloat& operator +=(const refloat& b);
    refloat& operator -=(const refloat& b);
    refloat& operator *=(const refloat& b);
    refloat& operator /=(const refloat& b);
    refloat operator -() const;
    
    __mpfr_struct* value;
  private:
    bool is_temp;
    bool is_set;
};

typedef std::vector<refloat> refloat_vec;


extern thread_local refloat *refloat_cache1;
extern thread_local bool refloat_all_allocations_temporary;
extern thread_local std::vector<__mpfr_struct*>* temporaries;


inline void setAllocationsToTemporary(bool tt){refloat_all_allocations_temporary=tt;}


void init_refloat_cache(mpfr_prec_t p);
void set_refloat_precision(mpfr_prec_t p);
void clear_refloat_cache();

void init_refloat_cache(void *arg);
void set_refloat_precision(void *arg);
void clear_refloat_cache(void *arg);

void refloat_finishUp(void *arg=NULL);
void refloat_startUp(void* arg=NULL);

void allocate_mpfr_struct(__mpfr_struct*& thestruct);
void init_mpfr_value(__mpfr_struct*& thestruct, mpfr_prec_t prec);

void swap(refloat& a, refloat& b);
void neg(const refloat& a, refloat& c);
void add(const refloat& a, const refloat& b, refloat& c);
void sub(const refloat& a, const refloat& b, refloat& c);
void mul(const refloat& a, const refloat& b, refloat& c);
void sqr(const refloat& a, refloat& c);
void div(const refloat& a, const refloat& b, refloat& c);
void sqrt(const refloat& a, refloat& b);
void abs(const refloat& a, refloat& b);
void mul_power_2(const refloat& a, long int p, refloat& c);
void sin_cos(const refloat& z, refloat& sn, refloat& cs);
void cos(const refloat& z, refloat& cs);
void sin(const refloat& z, refloat& sn);
void exp(const refloat& z, refloat& e);
void log(const refloat& z, refloat& l);
void atan2(const refloat& y, const refloat& x, refloat& v);
void erf(const refloat& z, refloat& e);
void pow(unsigned long int n, const refloat& x, refloat& v);
int compare(const refloat&a, const refloat& b);
int compare_abs(const refloat&a, const refloat& b);
void round(const refloat& a, refloat& c);
void Si_Ci_Acton(const refloat& z, refloat& Si, refloat& Ci, unsigned long int n);
void Si_Ci_Acton(const refloat& z, refloat& Si, refloat& Ci);
void Si_Ci_Continued_Fraction(const refloat& z, refloat& Si, refloat& Ci);
void Si_Ci_power_series(const refloat& z, refloat& Si, refloat& Ci);
void Si_Ci(const refloat& z, refloat& Si, refloat& Ci);
void set_pi(refloat& val);


inline void refloat::prefetch() const
{
  unsigned long int word_size = 8*sizeof(unsigned long int);
  unsigned long int nwords = mpfr_get_prec(value)/word_size;
  for(unsigned long int i=0;i<nwords;++i)
  {
    __builtin_prefetch(&((value->_mpfr_d)[i]));
  }
}


inline void allocate_mpfr_struct(__mpfr_struct*& thestruct)
{
  void* val = (void*)thestruct;
  if(posix_memalign(&val, 64, sizeof(__mpfr_struct)) != 0)
  {
    throw std::runtime_error("allocate_mpfr_struct : unable to allocate aligned memory");
  }
  thestruct = (__mpfr_struct*)val;
}


inline void init_mpfr_value(__mpfr_struct*& thestruct, mpfr_prec_t prec)
{
  if(prec < 16)
  {
    prec = mpfr_get_prec(refloat_cache1->value);
  }
  
  void* val = (void*)(thestruct->_mpfr_d);
  
  unsigned long int size = (((long int)(prec))>>6);
  size += 1;
  if(prec%64 != 0){size += 1;}
  size *= 8;
  
  if(posix_memalign(&val, 64, size  ) != 0)
  {
    throw std::runtime_error("init_mpfr_value : unable to allocate aligned memory");
  }
  (thestruct->_mpfr_d) = (mp_limb_t*)val;
  (thestruct->_mpfr_prec) = prec;
}


inline refloat::refloat(std::nullptr_t ptr) : is_temp(refloat_all_allocations_temporary), is_set(false)
{
  mpfr_prec_t prec = mpfr_get_prec(refloat_cache1->value);
  if(is_temp==false)
  {
    allocate_mpfr_struct(value);
    init_mpfr_value(value, prec);
  }
  else
  {
    set_temporary(prec);
  }
}


inline refloat::refloat(long int re, mpfr_prec_t prec, bool temporary) : is_temp(temporary || refloat_all_allocations_temporary), is_set(true)
{
  if( (prec < 16) && (mpfr_get_prec(refloat_cache1->value) > 16) ){prec = mpfr_get_prec(refloat_cache1->value);}
  if(is_temp==false)
  {
    allocate_mpfr_struct(value);
    init_mpfr_value(value, prec);
  }
  else
  {
    set_temporary(prec);
  }
  mpfr_set_si(value, re, GMP_RNDN);
}


inline refloat::refloat(const refloat& z) : is_temp(z.isTemp() || refloat_all_allocations_temporary), is_set(z.isSet())
{
  if(is_temp==false)
  {
    allocate_mpfr_struct(value);
    init_mpfr_value(value, mpfr_get_prec(z.value));
  }
  else
  {
    set_temporary(mpfr_get_prec(z.value));
  }
  if(is_set==true){set(z);}
}


inline refloat::refloat(std::string re, mpfr_prec_t prec, int base, bool temporary) : is_temp(temporary || refloat_all_allocations_temporary), is_set(true)
{
  if( (prec < 16) && (mpfr_get_prec(refloat_cache1->value) > 16) ){prec = mpfr_get_prec(refloat_cache1->value);}
  if(is_temp==false)
  {
    allocate_mpfr_struct(value);
    init_mpfr_value(value, prec);
  }
  else
  {
    set_temporary(prec);
  }
  if(re != ""){mpfr_set_str(value, re.c_str(), base, GMP_RNDN);}
}


inline unsigned long int refloat::to_ulong()
{
  return mpfr_get_ui(value, GMP_RNDN);
}


inline bool refloat::isTemp() const
{
  return is_temp;
}

inline bool refloat::isSet() const
{
  return is_set;
}


inline bool refloat::is_nan()
{
  return (mpfr_nan_p(value) || mpfr_inf_p(value));
}


inline void refloat::set(std::string re, mpfr_prec_t prec, int base)
{
  is_set = true;
  if((is_temp==false) || (!temporaries)){free(value->_mpfr_d);free( value );}
  else
  {
    temporaries->push_back(value);
  }
  if((is_temp==false) || (!temporaries))
  {
    allocate_mpfr_struct(value);
    init_mpfr_value(value, prec);
  }
  else
  {
    set_temporary(prec);
  }
  mpfr_set_str(value, re.c_str(), base, GMP_RNDN);
}


inline void refloat::set(std::string re, int base)
{
  is_set = true;
  mpfr_set_str(value, re.c_str(), base, GMP_RNDN);
}


inline void refloat::set(long int re, mpfr_prec_t prec)
{
  is_set = true;
  if((is_temp==false) || (!temporaries)){free(value->_mpfr_d);free( value );}
  else
  {
    temporaries->push_back(value);
  }
  if(is_temp==false)
  {
    allocate_mpfr_struct(value);
    init_mpfr_value(value, prec);
  }
  else
  {
    set_temporary(prec);
  }
  mpfr_set_si(value, re, GMP_RNDN);
}


inline void refloat::set(long int re)
{
  is_set = true;
  mpfr_set_si(value, re, GMP_RNDN);
}


inline void refloat::set(const refloat& x)
{
  is_set = true;
  mpfr_set(value, x.value, GMP_RNDN);
}


inline void refloat::set_prec(mpfr_prec_t prec)
{
  if(prec != mpfr_get_prec(value))
  {
    mpfr_set((*refloat_cache1).value, value, GMP_RNDN);
    free(value->_mpfr_d);free( value );
    allocate_mpfr_struct(value);
    init_mpfr_value(value, prec);
    mpfr_set(value, (*refloat_cache1).value, GMP_RNDN);
  }
}


inline mpfr_prec_t refloat::get_prec() const
{
  return mpfr_get_prec(value);
}


inline void swap(refloat& a, refloat& b)
{
  (*refloat_cache1).set(a);
  a.set(b);
  b.set(*refloat_cache1);
}


inline void refloat::neg()
{
  mpfr_neg(value, value, GMP_RNDN);
}


inline void neg(const refloat& a, refloat& c)
{
  mpfr_neg(c.value, a.value, GMP_RNDN);
}


inline void refloat::add(const refloat& x)
{
  mpfr_add(value, value, x.value, GMP_RNDN);
}


inline void refloat::add(long int i)
{
  mpfr_add_si(value, value, i, GMP_RNDN);
}


inline void add(const refloat& a, const refloat& b, refloat& c)
{
  c.set(a);
  c.add(b);
}


inline void refloat::sub(const refloat& x)
{
  mpfr_sub(value, value, x.value, GMP_RNDN);
}


inline void refloat::sub(long int i)
{
  mpfr_sub_si(value, value, i, GMP_RNDN);
}


inline void sub(const refloat& a, const refloat& b, refloat& c)
{
  c.set(a);
  c.sub(b);
}


inline void refloat::mul(const refloat& x)
{
  mpfr_mul(value, value, x.value, GMP_RNDN);
}


inline void mul(const refloat& a, const refloat& b, refloat& c)
{
  c.set(a);
  c.mul(b);
}


inline void refloat::sqr()
{
  mpfr_sqr(value, value, GMP_RNDN);
}


inline void sqr(const refloat& a, refloat& c)
{
  c.set(a);
  c.sqr();
}


inline void refloat::div(const refloat& x)
{
  mpfr_div(value, value, x.value, GMP_RNDN);
}


inline void div(const refloat& a, const refloat& b, refloat& c)
{
  c.set(a);
  c.div(b);
}


inline void refloat::sqrt()
{
  mpfr_sqrt(value, value, GMP_RNDN);
}


inline void sqrt(const refloat& a, refloat& b)
{
  b.set(a);
  b.sqrt();
}


inline void refloat::abs()
{
  mpfr_abs(value, value, GMP_RNDN);
}


inline void abs(const refloat& a, refloat& b)
{
  b.set(a);
  b.abs();
}


inline void refloat::mul_power_2(long int p)
{
  if(p>0)
  {
    mpfr_mul_2si(value, value, p, GMP_RNDN);
  }
  else
  {
    mpfr_div_2si(value, value, -p, GMP_RNDN);
  }
}


inline void mul_power_2(const refloat& a, long int p, refloat& c)
{
  c.set(a);
  c.mul_power_2(p);
}


inline void refloat::reciprocal()
{
  (*refloat_cache1).set((*this));
  this->set(1);
  this->div((*refloat_cache1));
}


inline void refloat::round()
{
  mpfr_round(value, value);
}


inline void round(const refloat& a, refloat& c)
{
  c.set(a);
  c.round();
}


inline void sin_cos(const refloat& z, refloat& sn, refloat& cs)
{
  mpfr_sin_cos(sn.value, cs.value, z.value, GMP_RNDN);
}


inline void cos(const refloat& z, refloat& cs)
{
  mpfr_cos(cs.value, z.value, GMP_RNDN);
}


inline void sin(const refloat& z, refloat& sn)
{
  mpfr_sin(sn.value, z.value, GMP_RNDN);
}


inline void exp(const refloat& z, refloat& e)
{
  mpfr_exp(e.value, z.value, GMP_RNDN);
}


inline void log(const refloat& z, refloat& l)
{
  mpfr_log(l.value, z.value, GMP_RNDN);
}


inline void atan2(const refloat& y, const refloat& x, refloat& v)
{
  mpfr_atan2(v.value, y.value, x.value, GMP_RNDN);
}


inline void erf(const refloat& z, refloat& e)
{
  mpfr_erf(e.value, z.value, GMP_RNDN);
}


inline void gamma(const refloat& z, refloat& e)
{
  mpfr_gamma(e.value, z.value, GMP_RNDN);
}


inline void pow(unsigned long int n, const refloat& x, refloat& v)
{
  v.set(1, mpfr_get_prec(x.value));
  (*refloat_cache1).set(x);
  
  while(n>0)
  {
    if(n&1)
    {
      v.mul((*refloat_cache1));
    }
    (*refloat_cache1).sqr();
    n = (n>>1);
  }
}


inline void set_pi(refloat& val)
{
  mpfr_const_pi(val.value, GMP_RNDN);
}


inline int compare(const refloat&a, const refloat& b)
{
  return mpfr_cmp(a.value, b.value);
}


inline int compare_abs(const refloat& a, const refloat& b)
{
  return mpfr_cmpabs(a.value, b.value);
}

#include "refloat_operators.h"



inline const refloat& conj(const refloat& x)
{
  return x;
}


inline const refloat& real(const refloat& x)
{
  return x;
}


inline refloat imag(const refloat& x)
{
  refloat temp(0, mpfr_get_prec(x.value), true);
  return temp; 
}


inline refloat abs(const refloat&  x)
{
  refloat temp(0, mpfr_get_prec(x.value), true);
  temp.set(x);
  temp.abs();
  return temp; 
}


inline refloat fabs(const refloat&  x)
{
  refloat temp(0, mpfr_get_prec(x.value), true);
  temp.set(x);
  temp.abs();
  return temp; 
}


inline refloat abs2(const refloat& x)
{
  refloat temp(0, mpfr_get_prec(x.value), true);
  temp.set(x);
  temp.sqr();
  return temp; 
}


inline refloat sqrt(const refloat& x)
{ 
  refloat temp(0, mpfr_get_prec(x.value), true);
  temp.set(x);
  temp.sqrt();
  return temp; 
}


inline refloat exp(const refloat&  x)
{ 
  refloat temp(0, mpfr_get_prec(x.value), true);
  exp(x, temp);
  return temp; 
}


inline refloat log(const refloat&  x)
{ 
  refloat temp(0, mpfr_get_prec(x.value), true);
  log(x, temp);
  return temp; 
}


inline refloat gamma(const refloat&  x)
{ 
  refloat temp(0, mpfr_get_prec(x.value), true);
  gamma(x, temp);
  return temp; 
}


inline refloat sin(const refloat&  x)
{ 
  refloat temp(0, mpfr_get_prec(x.value), true);
  sin(x, temp);
  return temp; 
}


inline refloat cos(const refloat&  x)
{ 
  refloat temp(0, mpfr_get_prec(x.value), true);
  cos(x, temp);
  return temp; 
}


inline refloat pow(const refloat& x, refloat y)
{ 
  refloat temp(0, mpfr_get_prec(x.value), true);
  log(x, temp);
  y.mul(temp);
  exp(y, temp);
  return temp; 
}


inline refloat pow(const refloat& x, unsigned long int n)
{
  refloat temp(0, mpfr_get_prec(x.value), true);
  pow(n, x, temp);
  return temp;
}


inline refloat pow(unsigned long int n, const refloat& x)
{
  refloat temp(0, mpfr_get_prec(x.value), true);
  pow(n,x,temp);
  return temp;
}














#endif
